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Abstract. I study the properties of the equihbrium probabihty distribution of a 
protein folding model originally introduced by Wako and Saito, and later reconsidered 
by Munoz and Eaton. The model is a one-dimensional model with binary variables 
and many-body, long-range interactions, which has been solved exactly through a 
mapping to a two-dimensional model of binary variables with local constraints. Here 
I show that the equilibrium probability of this two-dimensional model factors into 
the product of local cluster probabilities, each raised to a suitable exponent. The 
clusters involved are single sites, nearest-neighbour pairs and square plaquettes, and 
the exponents are the coefficients of the entropy expansion of the cluster variation 
method. As a consequence, the cluster variation method is exact for this model. 



1. Introduction 

The cluster variation method (CVM) is an approximate method of equihbrium statistical 
physics introduced by Kikuchi pOj as a generalization of the Bethe-Peierls |21 E] and 
Kramers-Wannier jH E] approximations. In its modern formulation jS] it is based on 
the minimization of an approximate variational free energy which is derived from the 
exact one by a truncation of the cumulant expansion of the entropy. This free energy 
depends on the probabihty distributions of local clusters. An account of the CVM and 
its applications up to the beginning of the '90s can be found in 0. A review article is 
also in preparation by the present author 

Recently the relationship between the Bethe-Peierls approximation, that is the 
lowest order CVM approximation, and the belief propagation method P, widely used 
for inference and optimization problems defined in terms of probabilistic graphical 
models, has been discovered [101 IH] , and this has led to the development of the so- 
called generalized belief propagation (GBP) jTO], an iterative message-passing algorithm 
whose fixed points are stationary points of the CVM free energy. 

There are only a few cases in which the CVM, and hence the GBP algorithm, gives 
exact results. In most cases the exactness is due to the tree-like topological structure 

§ alessandro.pelizzola@polito.it 



Exactness of CVM for a protein folding model 



2 



of the underlying lattice or graph. Leaving apart tree-like structures, the only case for 
which I am aware of the exactness of the CVM is that of disorder varieties ^21 IT^ . 
which occur in two-dimensional Ising models with short-range competitive interactions 
on ordinary translation-invariant lattices, where the correlations have a particularly 
simple form which resembles that of one-dimensional models. 

In all the cases in which the CVM is exact, independent of the reason, the 
equilibrium probability distribution factors into a product of local cluster probabilities, 
each raised to a suitable integer exponent. 

Given the limited number of cases in which the CVM gives an exact solution, it 
is particularly relevant to show that for a particular protein folding model the CVM is 
exact and the equilibrium probabihty factors. This is the purpose of the present paper. 

The model I am going to study has been introduced and exactly solved for the first 
time in 1978 by Wako and Saito [Ellin], who showed its potential apphcation to the 
protein folding problem. Much later, Munoz and Eaton El E| reconsidered this 
model, in a slightly different formulation, solved it within the so-called single (double, 
triple) sequence approximation, and made comparisons with their own experimental 
data, after which the model became rather popular (see for instance PH])- Flammini, 
Banavar and Maritan [20] derived a mean-field theory, made Monte Carlo simulations 
and solved exactly the a-helix case and later Bruscolini and I [211 122] solved exactly 
the general case, using ideas that are essentially the same as those by Wako and Saito 
|14[ [T3] . The same model and techniques have recently been used for strained epitaxy 
on a modulated substrate [^ . 

In all these papers no reference was ever made to the original works by Wako and 
Saito. As far as I know, the first paper citing both Wako-Saito and Munoz-Eaton is 
due to Itoh and Sasai [23, which refer to the model as the Wako-Saito-Munoz-Eaton 
(WSME) model, as I'll do in the following. 

The WSME model is a one-dimensional effective model, where a given protein is 
regarded as a sequence of monomers (residues) connected by peptidic bonds, and a 
binary variable is associated to each peptidic bond, with the value 1 representing the 
folded (or native) configuration and the value representing the unfolded configuration. 
An entropic cost is associated to the folded configuration. Only native interactions are 
considered, like in Go [^ models, that is two residues are allowed to interact only if 
they are in contact in the native state. Moreover, and this is the main characteristic of 
the present model, two residues can interact only if all the peptidic bonds between them 
in the sequence are in the native state, which gives rise to the many-body, long-range 
interactions. 

The exact solution of the equilibrium thermodynamics of the model can be obtained 
through a mapping to a two-dimensional model [25 , defined on a triangular-shaped 
portion of the square lattice, with local constraints as the only interactions. The 
dimension of the state space of a row of this model is at most equal to the length 
of the protein and this makes the transfer matrix approach feasible. 

In the present paper I show that the equilibrium probabihty distribution of this 
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model in its two-dimensional version can be written in a factorized form which is 
typical of generalized mean-field theories. It is a product of local cluster probabilities, 
marginals of the global probability. In this product only single-site, nearest-neighbour 
(NN) and square plaquette clusters appear, and the corresponding probabilities are 
raised to exponents which are the coefficients of the entropy expansion which appears 
in the formulation of the CVM 12]. As a consequence, the CVM turns out to be exact 
(which was already noticed empirically in |2I]), providing an exact variational free energy 
depending on a number of variables which scales only quadratically with the length of 
the protein. 

The plan of the paper is as follows: in Section El the model is described in detail, its 
mapping to a two-dimensional model is discussed in Sectional then the factorization of 
the probability distribution is shown in Section HI The CVM is described in Section 
where its exactness for the present model is shown and the consequences of this property 
are discussed. Finally, conclusions are drawn in Section IHl 

2. The model 

The WSME model, if considered as a model for protein folding, is a simplified effective 
model in the sense that one considers a very restricted state space, and imagines to 
integrate with respect to all the degrees of freedom except the binary peptidic bond 
variables, obtaining a Hamiltonian which is actually an effective free energy. 

In order to describe the model, consider a protein of length + 1 residues, and 
associate a binary variable mj, i = 1,. . .N to each peptidic bond between residues i 
and i + 1. = 1 denotes a peptidic bond in a native configuration, mj = a peptidic 
bond in an unfolded one. Since there are actually many more unfolded configurations 
than native ones, an entropic cost Asj < is associated to each residue. 

Two residues are allowed to interact only if they are in contact in the native state. 
A detailed definition of contact between residues is not needed here and can be found 
in [TBI [T71 ITH] . Here I simply assume that a contact matrix A is given, with elements 
Aij = 1 if residues i and j + I (or, equivalently, peptidic bonds i and j) are in contact 
in the native state, and otherwise. 

As a further condition, residues i and j + 1 interact (with an energy eij) only if all 



where T is the temperature. 

The remainder of this paper does not focus on the applications to the protein folding 
problem, but rather on the mathematical properties of the statistical mechanics of the 
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Figure 1. A typical configuration of the WSME model. An empty (resp. filled) circle 
at row j, column i represents the variable rtij^i taking value (resp. 1). 

model itself, therefore from now on we shall deal with a generic Hamiltonian of the form 

N N j 

H{m)=Y.Y.h^,Y[^k, (2) 

i=l j=i k=i 

where the hi/s can be temperature dependent. 
3. Mapping to a 2— dimensional model 

j 

The above Hamiltonian suggests to introduce the new binary variables rnj^i = J^m^ 

k=i 

(notice the order of the indices), which take value 1 if = 1, i < k < j , and otherwise. 
These new variables can be associated to the triangular-shaped portion of the square 
lattice defined by 1 < i < j < and shown in Figure ^ and the original variables are 
included in this set, since j = mj. 

The Hamiltonian can now be written as 

N j 

H{m) = J2Y1 Kj'^j.h (3) 
j=i i=i 

but the new variables are not all independent. Since we have a model with 2^ 
configurations written in terms of N{N + l)/2 binary variables, constraints must exist 
between these variables. These constraints can be written in different ways, and we 
chose 121] to write them as 

m^. . = mj_i,imj-j+i, I <i < j < N. (4) 
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One is therefore left with a model with a local (actually non-interacting) 
Hamiltonian, and local constraints in place of the interactions. Formally one could 
also include the constraints into the Hamiltonian, by studying the limit A — +00 of 

N j N j-1 

H{m) = X] + ^Y1 ~" (5) 

i=l 1=1 j=2 i=l 

In the new representation the feasibility of the transfer matrix approach becomes 
evident. Consider the variables in row j in Figure ^ that is rrij = {mj^i,mj^2, ■ ■ -f^jj)- 
Because of the constraints all configurations with mj i = 1 and mj j+i = are forbidden, 
leaving as the only allowed states for rrij the j + 1 states denoted by and defined by 

(0, i<k 

, < A; < J. (6) 

In Figure HI these states are readily depicted by placing k O's on the left of row j and 
j — k I's on the right. 

A transfer matrix solution for a protein of length + 1 has then to deal with the 
product of — 1 matrices, the largest of which has size (A^ + 1) x A^, the number of 
operations involved grows polynomially in A^, and the model can be easily solved even 
for long proteins. 



4. Factorization of the probability distribution 

The existence of a row-to-row transfer matrix is intimately connected to the 
factorization of the model's probability distribution over rows. This is well known in 
statistical mechanics (see for instance |2E])- Readers interested in a purely statistical 
formulation will regard this as an instance of the junction tree theorem j2Zll2H]- I report 
here the basic steps of a proof for the present case. 

First of all, write the Hamiltonian as a sum of row-to-row terms: 

7V-1 

Him) = J2 Hjj+i{mj,mj+i). (7) 
i=i 

If the Hamiltonian has the form reported in Equation (jSJ, then one can write 

j i+i j 

Hjj+i = bj hijnij^i + bj+i ^ hij+iruj+i^i + A ^(m^+i^i - mj^imj+i^i+iY , (8) 

1=1 i=l 1=1 

where, in order to take into account the boundaries, bi = = 1 and bj = 1/2 for 
1 < j < A^. For simplicity of notation, I drop from now on the arguments of Hj j^i. 
Of course, the specific form of -f^jj+i is irrelevant here, only the splitting into terms 
involving only adjacent rows is important. 

Then we introduce the Boltzmann distribution 
1 1 

p{m) = — exp(-i^) = — exp(-ifi,2 Hn_i^n), (9) 



Exactness of CVM for a protein folding model 



6 



where we have absorbed f3 = (kBT) ^ in the definition of the Hamiltonian and 
introduced the partition function 

Z = J2exp{-H) = J2exp{-H,,2 Hn-i,n). (10) 

m m 

Finally we define the marginal probability distributions for a row 

pA^j)= E pM' ^<j<n, (11) 

{mi,lj^j} 

and for a pair of adjacent rows 

Pi,i+i(^i> ^i+i) = E P("^)' 1 < j < - 1- (12) 

{■mi,lytj,j+l} 

The arguments of these probability will also be dropped in the following. 

Introducing the partial partition functions for upper and lower portions of the 
lattice 

Zujiruj) = ^M-Hj-ij) ■ ■ ■ Eexp(-ifi,2), 

rrij - 1 mi 

Zijimj) = J2 exp(-i/j, j+i) ■ ■ ■ ^ exp(-i7jv-i,^), 

mj^l m^ 

together with the boundary conditions 

ZuAmi) = 1, 

ZL,NimN) = 1, 
it is easy to show that 

Pj,j+i = ^Zuj{mj)ZLj+i{mj+i)exp{-Hjj+i), 

and hence 



(13) 



(14) 



(15) 



Pl,2- ■ ■PN-1,N ^ f Tj\ f \ nci\ 

= — exp[—H) = p[m). (16) 

P2---PN-1 Z 

This shows that the global probability can be written as a product of row-pair and row 
probabilities, raised to exponents 1 and -1 respectively, which is the first step towards 
our final result. As a consequence, a CVM with all row-pairs as maximal clusters (see 
Sectional) would be exact for this model, but this is not interesting here, since a much 
stronger property will be proved below. 

The above discussion is not specific to the present model and so far the constraints 
have not yet been exploited. The next step is then to make use of the constraints to 
write in a factor form similar to Equation (fTBj) the row and row-pair probabilities. 

Consider first the row probability. Observe that, due to the constraints Equation 

©, 

(17) 

ruj^i = 1 =^ ''^j,k = 1 \fk > i, 
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and the probability of the state defined in Equation (jHI) reduces to a NN pair 
probabihty: 



p,(ef)=pf+^(0,l), l<k<j 



where fc, ^+1) is the NN (horizontal) pair probability for the pair at row j 

and columns k,k + l. Indeed, is the only state of row j with ^ = and mj^k+i = 1- 
In addition, for the states (all I's) and e^- (all O's) one can write 



P.(e;)=pf(l,l), 



(19) 



(actually the probabilities for these states could be reduced to site probabilities, but 
this is not useful here). 

Introducing the site probabilities p*(mj j) and observing that 

pr'(o,o)=pf^(o), 

one can write the above results for the row probability in the /c-independent factorized 
form 

^.^gfc) ^ pj''(m,,i, mj^2) ■ ■ -pf '''(m,-,,!, m,-,) ^^^^ 

A similar result can be obtained for the row-pair probability. In this case one 
needs to define the square probability pjj^\(mj^j, mj_j+i; mj+i^i, mj+i^j+i), the triangle 
probability p'j-'j1i^\{mjj] mj+ij, mj+ij+i) for the triangles lying on the diagonal boundary 
and the NN (vertical) pair probability p*j_,_]^(mj^j; mj+i^i). 

There are two kinds of row-pair configurations: {ej,ej^j) which represents the 
cases in which peptidic bonds j and j + 1 are in the same native stretch (e.g. rows 3 and 
4 in Figure or (for k = j) a new native stretch starts at j + 1 (e.g. rows 7 and 8 in 
Figure QJ, and (e^, e^-^}) which represents the cases in which a native stretch of length 
j — k ends at j (e.g. rows 5 and 6 in Figure P). For both kinds we have to show that a 
factorization like Equation (j?H) occurs. 

The proof parallels the one for the row probability. Consider the first kind of 
configurations and observe that, due to the constraints, 

P,,M^^, e';^,) = pfj:l:l{o, i; o, i), i<k< j. (22) 

For the boundary cases one has 

Pjj+i(e°, e°+i) = pj;|+i(l, 1; 1, 1), 
Pj,j+iie], ej+i) = pfl^lliO; 0, 1). 
Observing that 

p;;^(0,0;0,0)=p}+^i(0;0) 
p;;^(l,l;l,l)=p},^+i(l;l) 



(23) 



(24) 
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the row-pair probability can be written as 

4+1) = ^^^'"^'^"f^^ 0<k<j, (25) 
Pj,j+i ' ' ' PjJ+i 

where the last factor in the numerator is a triangle probability. For simplicity, the 
arguments of the probability have been dropped. 

Consider now configurations of the second kind. Equations ()22j] . ()23|1 and ()25j) are 
replaced by 

p,-,+i(eJ, e^l) = p]^!;tl{^, 1; 0, 0), l<k< j, (26) 



p,,,+i(eO,e^:J)=pj;^i(l,l;0,0), 

p..+i(e^,e^:})=P^;Sl(o;o,o) 



(27) 



and 



1,2 

e^::}) = ^^'^^ ^^.^^^^^^^^^ < A; < (28) 
^'i ■ ■ ■ Pj,j+i 

respectively. 

The factorization property 

p,,^.i(m,, m,+0 = P^^-^-'-P^^^^jP^Sl o<k<j (29) 

■ ■ ■ Pj,j+i 

is then proved for any valid row-pair configuration. 

In order to obtain the final result of this section, plug Equations (PT|) and into 
Equation (fTH|) . obtaining 

^A^-l j \ (N-lj-l 

nm:S\ nn^} 

P{m) = -VlZiiz^ Jr'';' \ , (30) 

n m^+i n m^' 

j=2 i=2 I \j=2 1=1 

that is the global probability is reduced to a product of square, triangle and site 
probabilities divided by a product of pair probabihties. 



5. Exactness of the cluster variation method 



In this section, after a brief introduction to the CVM, it is shown that the factorization 
Equation (jHUjl implies the exactness of the CVM, and the consequences of this property 
are discussed. 

The CVM in its modern formulation is derived from the variational principle 
of statistical mechanics, which states that, given a model with variables m and 
Hamiltonian H{m), its equilibrium (Boltzmann) probability distribution peq(^) = 
exp(— i7(m))/Z is the distribution which minimizes the variational free energy 

F = U — S = '^p{m)H{m) + ^p(rn,) Inp(m). (31) 

m m 
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The physical values of the free energy F, the energy U and the entropy S can then be 
calculated at the minimum. 

If the variables are associated to the nodes of a graph, or to the sites of a lattice, 
one can define clusters of nodes (e.g. single sites, pairs, triangles, square plaquettes, 
. . . ). For each cluster a define also the corresponding set of variables rUa, probability 
distribution Pa{ma) and entropy 

S'a = -^]9„(ma)lnp„(m„). (32) 
The entropy cumulants are defined by 

Sa=J2 SfS, (33) 
f3Ca 

which can be solved with respect to the cumulants by means of a Mobius inversion, 
which yields 

S(,= J2i-ir--^^S^, (34) 

aC/3 

where Ua denotes the number of nodes in cluster a. 
The full entropy 

S = Y.Sp (35) 

can then be approximated by selecting a set R of clusters, made of certain maximal 
clusters and all their subclusters, and truncating the cumulant expansion by retaining 
only terms which correspond to clusters in R. One obtains 

S-Y.Sf3=Y, aaSa, (36) 

where the coefficients a^, sometimes called Mobius numbers, satisfy |H] 

aa = l ypeR. (37) 

The above condition, practically useful for determining the a^s, means that every 
subcluster must be counted exactly once in the entropy expansion. 

Now assume that the Hamiltonian is made of local terms only, and R has been 
chosen such that one can write 

H{m) = J2 Ha{m^). (38) 

aeR 

In such a case the variational free energy can be written, with the above approximation 
on the entropy and no approximation on the energy, as 

FcvM = YPo^('^o)H^{mo) + Y ciaYP^'^^) lnp(m,c,), (39) 

oS-R m.a a£R rria 

where the minimization must be performed with respect to the p^'s with the constraints 

YPai'TTla) = 1, Vft G R, 

ma {4:0) 

Y Pc{rna) = Pfsimp), V/3 C « G i?. 

ma\i3 
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It is interesting to observe that a suitable factorization of the equihbrium probabihty 
of the model implies that the variational free energy Equation ()31|) reduces to the CVM 
variational free energy Equation ()H9|1 with no approximation. More precisely, assume 
that the equilibrium probability factorizes in terms of its marginals according to 

Peq(m) = n [Pa{m^)r ■ (41) 

Qg-R 

Then one can restrict the variational principle to distributions p{m) with the same 
property. The entropy of such a distribution is 

S = — ^]9(m) Inp(m) = 

m 

= -J2aaJ2p{m) In pairria) = (42) 

a m 
a rricy, 

and the CVM free energy is therefore obtained with no approximation (recall that no 
approximation was made on the energy term). 

The WSME model in its two-dimensional representation falls precisely in this case. 
If one chooses as maximal clusters all the square plaquettes and the triangles lying on 
the diagonal boundary the Hamiltonian can obviously be written as in Equation ()38p 
and the Mobius numbers for the entropy expansion, obtained by Equation pTjl . are: 

• Oq, = 1 for the maximal clusters; 

• aa = for the triangles not lying on the diagonal boundary and the next -nearest- 

neighbour pairs, since they are contained in exactly 1 maximal cluster; 

• eta = — 1 for all the NN pairs contained in 2 maximal clusters, that is all NN pairs 

except the boundary ones; 

• Cq, = for the boundary NN pairs; 

• = 1 for the sites contained in 4 maximal clusters and 4 NN pairs, that is all sites 
except the boundary ones; 

• Oq, = for the boundary sites. 

With the above Mobius numbers it is immediate to check that the factorization 
Equation (PT|) is exactly the one that we have obtained in Equation (jHUI) for the WSME 
model, and hence the CVM is exact for this model. The corresponding variational free 
energy, dropping all the arguments in the entropy term, reads 

N j N-1 j 

j=l i=l rrij^i j=l i=l 

N~l j N-lj-1 

- E ep}.+i - E Y.pr ^pr+ (43) 

j=2 1=2 j=2 1=1 

N-lj-1 

+ E EpJi^Pi- 

j=3 i=2 
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The local constraints Equation (@]) can be either included in the Hamiltonian or imposed 
by hand on the probabilities. 

The numerical minimization of the above variational free energy can in principle be 
performed by a provably-convergent, double-loop algorithm like the one proposed by 
Heskes, Albers and Kappen 12^, while the GBP fails, probably due to the constraints, 
except in very simple cases. This is not an important point however, since the strength 
of this result is not that it improves on the transfer matrix method, which is already 
very efficient in solving for the equilibrium. On this side the only advantage of this 
approach is probably that in this scheme it is easier to calculate correlation functions, 
since the probabilities are directly accessible. 

The real strength of this result, apart from having found a new model which is 
exactly solvable by the CVM, is that here the CVM can serve as a basis to build 
powerful approximation for the dynamics of the model, which is extremely relevant for 
the protein folding problem. For this purpose, it might also be useful to observe that 
the variational free energy can be written explicitly in terms of the local expectations 

Xj,i = {ruj^i) = ^j,iP]i^j,i) =P]i^)- (44) 

1=0,1 

To this end observe that in the variational free energy Equation ()43|) the energy term is 
already a linear function of these expectations, while the probabilities appearing in the 
entropy term can be written as functions of these expectations as independent variables. 
For instance, using Equation and normalization one obtains 

For the NN pair one configuration is forbidden by the constraints, and the 
remaining probabilities are determined by normalization and marginalization to the 

site probabilities, obtaining 

1,1+1 



pf (0, 0) = 1 - Xj^i+i 

pf^^{^, 1) = XjA+i - Xj^i (46) 



p'-'+^l, 1) = 



for the horizontal pair and 

p;j+i(o,o) = i-Xj^i 

0) = Xj^i - Xj+i^i (47) 

for the vertical pair. Similarly, for a triangle lying on the diagonal boundary, only 
4 configurations are allowed by the constraints, and their probabilities are again 
determined by normalization and marginalization to subclusters, with the result 



p}'j^\(0; 0, 0) = 1 - Xj.i - Xj+i^i+i + Xj+i^i 



0; l) ~ Xj+i^i+i — Xj+i 

Pj^^j^iil] 0, 0) = Xj^i — Xj^i^i 
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Finally, for a square plaquette we have 5 allowed configurations and the probabilities 
p5;^(0,0;0,0) = 1-x,,+i 
Pjj-+i(0, 1; 0, 0) = Xj^i^i + Xj^i^i — Xj^i — 

P^j+ii^, 1; 0, 1) = - (49) 

1; 0, 0) = Xj^i — Xj+i^i 

Substituting the above probabilities into the variational free energy Equation (PHj) yields 
an exact variational free energy as a function of N{N + 1)/2 independent local variables. 

6. Conclusions 

I have shown that the equilibrium probability factors into a product of local cluster 
probabilities, and hence the CVM is exact, for the WSME model of protein folding. After 
the one-dimensional WSME model, which has long-range, many-body interactions, has 
been mapped into a two-dimensional model which has local constraints as the only 
interactions, the proof goes through two steps. The first step exploits the locality of the 
interactions, the second one the detailed form of the constraints. 

The result is especially relevant on the methodological side, since leaving apart 
tree-like models, the CVM is exact, as far as I know, only on disorder varieties of 
two-dimensional models. 

Moreover, some consequences can be expected also on the model side. As far as 
equilibrium properties are concerned, we have almost no improvement with respect 
to the transfer matrix method, which is already very fast in this matter. The main 
advantage of the CVM solution is that correlation functions are much easier to calculate. 
The most important point is however that the CVM solution for the equilibrium can be a 
starting point for good approximations for the dynamics, which is of utmost importance 
in the context of the protein folding problem. Two approaches can be followed [SOI EI]- 
On one hand, one can develop a master equation approach based on the approximation 
that the state of the system at any time can be described as an equilibrium state of 
the WSME model, with a time-dependent Hamiltonian (the so-called local equilibrium 
approximation jSZ|)- On the basis of the results reported here, this means that the 
probability at any time is assumed to factorize as the equilibrium one. Alternatively, 
one can build an approximation for the dynamics with the path probability method 
(PPM, the dynamical version of the CVM) [Zj, with square plaquettes and triangles 
as maximal clusters. Since the stationary state of the PPM corresponds to the CVM 
solution for the equilibrium, one obtains an approximation for the dynamics which is 
guaranteed to converge to the exact equilibrium state. It can be verified jSOJ EI] that 
the two approaches are equivalent in the limit of vanishing time step and that very good 
agreement is obtained with respect to the exact solution for short chains. 
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